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Thermal  management  of  a  solid  oxide  fuel  cell  (SOFC)  stack  essentially  involves  control  of  the  tempera¬ 
ture  within  a  specific  range  in  order  to  maintain  good  performance  of  the  stack.  In  this  paper,  a  nonlinear 
temperature  predictive  control  algorithm  based  on  an  improved  Takagi-Sugeon  (T-S)  fuzzy  model  is  pre¬ 
sented.  The  improved  T-S  fuzzy  model  can  be  identified  by  the  training  data  and  becomes  a  predictive 
model.  The  branch-and-bound  method  and  the  greedy  algorithm  are  employed  to  set  a  discrete  optimiza¬ 
tion  and  an  initial  upper  boundary,  respectively.  Simulation  results  show  the  advantages  of  the  model 
predictive  control  (MPC)  based  on  the  identified  and  improved  T-S  fuzzy  model  for  an  SOFC  stack. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  solid  oxide  fuel  cell  is  a  high-  or  intermediate-temperature 
fuel  cell,  which  operates  in  the  range  600-1000  °C.  The  SOFC  system 
is  generally  considered  to  suit  the  generation  of  electricity  and  heat 
in  industrial  application.  A  high  operating  temperature  allows  inter¬ 
nal  reforming  and  improves  the  reaction  kinetics  [1  ].  However,  the 
fuel  cell  can  be  damaged  by  high  temperature  due  to  thermal  fatigue 
or  thermal  cracking.  High  temperature  also  affects  the  stack  relia¬ 
bility  and  durability,  and  shortens  the  stack  lifespan  [2].  Therefore, 
thermal  management  that  controls  the  operating  temperature  and 
reduces  temperature  fluctuation  is  very  important  for  SOFC  stacks. 

Many  offline  model-based  control  methods  have  been  estab¬ 
lished  for  fuel  cells  [3-9],  but  few  papers  deal  with  MPC  based  on 
online  control.  MPC  has  been  an  active  field  of  research  during  the 
last  three  decades,  and  there  have  been  numerous  successful  appli¬ 
cations  of  MPC  technology  [10-12].  MPC  [13,14]  has  a  number  of 
advantages  in  that  it  can  handle  multivariable  system  problems, 
allows  operation  states  close  to  the  constraints,  is  capable  of  aware¬ 
ness  of  the  actuator  limitations  of  the  model-based  control,  etc. 
MPC  consists  of  model  prediction,  receding  horizon  optimization, 
and  online  feedback  correction  [15].  The  SOFC  system  is  a  nonlinear 
one,  in  which  the  parameters  vary  within  an  operating  range,  and 
an  MPC  algorithm  can  satisfy  the  requirements  of  a  control  strategy. 
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Some  researchers  have  examined  control  methods  and  strate¬ 
gies  in  relation  to  the  SOFC  system.  Jurado  [16]  built  a  predictive 
control  model  to  achieve  online  control  of  an  SOFC  system,  which 
used  a  fuzzy  Hammerstein  model  identified  by  the  input-output 
data.  Wu  et  al.  [17]  applied  a  nonlinear  model  predictive  control 
method  in  order  to  control  the  voltage  and  guarantee  fuel  utiliza¬ 
tion  within  a  safe  range.  The  nonlinear  predictive  controller  was 
based  on  an  improved  radial  basis  function  (RBF)  neural  network, 
and  a  genetic  algorithm  (GA)  was  used  to  optimize  the  parameters. 
Yang  et  al.  [18]  presented  a  nonlinear  predictive  control  algorithm 
based  on  a  T-S  fuzzy  model  for  a  molten  carbonate  fuel  cell  stack. 

In  this  study,  an  online  nonlinear  MPC  scheme  based  on  an 
improved  T-S  fuzzy  model  has  been  built  to  control  the  temper¬ 
ature  within  a  safe  range.  This  improved  T-S  fuzzy  model  could 
be  identified  by  training  data  that  were  provided  by  a  physical 
model  from  Ref.  [19].  The  control  sequence  could  be  discretized 
and  optimization  could  be  sought  using  a  principle  of  the  branch- 
and-bound  method,  and  a  greedy  algorithm  has  been  employed 
to  set  the  initial  upper  boundary  for  the  performance  index  [20]. 
The  MPC  algorithm  met  the  requirements  for  an  online  tempera¬ 
ture  control  strategy,  and  the  simulation  results  reflected  the  more 
effective  temperature  control. 

2.  Structure  of  the  MPC  system  for  an  SOFC 

2.1  Temperature  predictive  control  system 

The  frame  of  the  nonlinear  MPC  system  of  an  SOFC  is  shown  in 
Fig.  1,  which  mainly  consists  of  a  controlled  plant  (SOFC  stack),  a 
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Fig.  1.  The  structure  of  the  MPC  system. 


nonlinear  improved  T-S  fuzzy  predictive  model,  and  a  predictive 
controller.  The  SOFC  physical  model  (see  Appendix  B)  described  in 
Ref.  [19]  was  built  with  MATLAB  to  imitate  the  real  3.5  kW  SOFC 
stack  and  is  shown  in  Fig.  2.  The  physical  model  was  constructed 
from  the  mass  balance  equations  and  the  enthalpy  balance  equa¬ 
tions.  In  addition,  the  MPC  system  also  contained  modules  of  the 
branch-and-bound  algorithm,  the  greedy  algorithm,  and  the  tapped 
delay  line  (TDL).  For  the  SOFC  dynamic  physical  model,  tempera¬ 
ture  was  the  controlled  variable.  Flow  rates  of  hydrogen  and  air 
were  chosen  as  manipulated  variables,  and  the  load  current  was 
considered  as  a  disturbance. 

In  Fig.  l,yre/  is  the  reference  curve  of  stack  temperature,  u(/<)  is 
the  manipulated  variable,  y(k)  is  the  current  value  of  stack  temper¬ 
ature,  and  y(k  +  1 )  is  the  predicted  next  value.  The  T-S  fuzzy  model 
predicts  future  stack  temperature  values  on  the  basis  of  history 
stack  temperature  values.  The  optimal  controller  can  determine 
the  next  control  signal  for  the  SOFC  stack  according  to  the  differ¬ 
ence  between  the  next  reference  values  and  the  predicted  stack 
temperature  values. 

2.2.  Predictive  model  based  on  an  improved  T-S  fuzzy  model 

An  improved  T-S  fuzzy  model  can  be  obtained  by  antecedent 
and  consequent  identification.  The  fuzzy  rule  is  of  the  form  “if. . . 


then. . .”.  The  ith  rule  of  the  Ith  output  yti{k  +  1 )  is  given  by 
Rt  i  :  If  x{k)  is  Al,  then 

9l,i(k  +  1)  —  P^o  +  Pj,i*kl  3 - (t  =  1,  .  .  •  ,  c)  (1) 

where  x{k)  is  the  regression  data  vector  consisting  of  input-output 
data  at  the  kth  instant  and  before,  c  is  the  number  of  rules,  Al  = 
Ki-  •••’<„>  is  the  set  of  membership  functions  associated  with 
the  ith  rule,  and  pi  =  [p!  0,  pi  1 , . . . ,  pi  ]  is  the  parameter  vector  of 
the  ith  submodel. 

Antecedent  identification  is  implemented  by  fuzzy  clustering 
based  on  the  principle  of  the  fuzzy  C-means  (FCM)  algorithm  [21  ]. 
The  consequent  part  of  the  fuzzy  rule  is  identified  by  using  the 
traditional  linear  identification  method,  the  least-squares  method. 

The  improved  T-S  fuzzy  model  is  fundamentally  different  from 
the  modified  T-S  fuzzy  model  in  Ref.  [19]  because  the  latter  is 
an  offline  model  method  whereas  the  former  is  an  online  model 
method.  Offline  models  are  sometimes  made  up  of  many  control 
rules  and  the  open-loop  data  that  they  generate  do  not  cover  all 
data  needed  for  closed-loop  control,  hence  it  is  difficult  to  obtain 
accurate  predictive  output  with  such  models. 

The  specific  algorithm  and  formula  of  the  improved  T-S  model 
is  inferred  as  follows  [21-25]: 

(1)  Initialization  of  the  number  of  fuzzy  clustering  c  (c  is  greater 
than  1),  the  clustering  center  vector  vf  =  [vn  . . .  vin]  (i  =  1,  . . ., 
c),  the  fuzzy  parameter  m  (m  is  greater  than  1),  the  discarding 
index  q(q  is  greater  than  1 ),  and  the  learning  rate  q{rj  is  greater 
than  0). 

(2)  Assessment  of  whether  there  is  a  clustering  center  vector  v *  in 
the  q  (discarding  index)  consecutive  sampling  period  that  sat¬ 
isfies  the  criterion  that  the  distance  from  to  the  consecutive 
input-output  data  vector  is  always  the  farthest.  If  so,  a  new  cen¬ 
ter  vector  is  adopted  to  replace  vif  for  example,  the  current  data 
vector  x(/<). 

(3 )  Calculation  of  the  distance  between  x(k)  and  each  cluster  center 
at  the  (k  -  1  )th  instant: 


n 

dW=>J2[X *(k)-¥k)]2.  '■  =  !>•..,  C  (2) 

N J=1 


Fig.  2.  The  SOFC  stack  physical  model. 
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and  accumulating  the  number  of  clustering  satisfying  the  crite¬ 
rion  that  the  distance  between  each  clustering  center  vector  and 
the  consecutive  input-output  data  vector  is  always  the  farthest. 

(4)  Determination  of  the  membership  degree  of  x(/<)  with  respect 
to  the  ith  clustering  at  the  (/<  -  1  )th  instant: 

-i  -i 


= 


(2/(m-l)) 


i  =  1 , . . . ,  c 


(3) 


(5)  Adjustment  of  the  clustering  center  vector  at  the  kth  instant 
according  to  the  membership  degree  ofx(/<)  with  respect  to  the 
ith  clustering  at  the  (/<  -  1  )th  instant: 


vi(k)  =  v,(k  - 1 )  ■ +  m'iikf  [x(k)  -  Vf(k  - 1 )] 


(4) 


(6)  Calculation  of  the  distance  between  x(/<)  and  each  clustering 
center  as  well  as  the  membership  degree  of  x(/<)  according  to 
the  new  clustering  center  at  the  kth  instant: 


d,(k)  = 


-  v#)]2 

\  J=1 


MiW  = 


E 

j=  i 


djWj 


(2/(m-l)) 


(5) 


(6) 


(7)  Use  of  the  least-squares  recurrence  method  to  identify  the 
linear  composition  factors  of  consequence.  The  data  vector 
x(fc)  =  [x^  xk2  . . .  Xfe„]  is  given.  The  output  of  the  fuzzy  model, 
y(k  +  1),  is  then  calculated: 


y(k  + 1)  = 


Mk  +  D 


c  c 

=  ^2^i(k)yi(k  +  1)  =  y>(feXp,o+P.xki  +  ■  ■  •  +  Pi„xkn) 

i=l  i=l 

=  <£(/<)  x  6> 

where  wI[x(/<)]  expresses  expectance  value, 

^lx(fc)] 


(7) 


Mi(fc)  = 


<ZW=[Mi  M2-  ■  Me MiXfci  M 2Xfci •  ■  McXU ■  •  -MiXto M2Xkn McXkn] 


0  =  [PlO  P20'  •  PcO Pll  P21 ' '  Pel  •  •  Pin  P2n-  •  Pen  ]T 

The  least-squares  recurrence  algorithm  is  employed  to  iden¬ 
tify  0  in  order  to  minimize  the  total  cumulative  square-root 
error  between  the  actual  output  and  the  model  output. 

(8)  If  the  control  process  does  not  end,  k  =  k-  1,  it  returns  to  Step  2. 

Completion  of  the  above  modeling  process  leads  to  a  nonlinear 
improved  T-S  fuzzy  predictive  model. 

3.  Discrete  optimization  of  control  variables 


In  this  study,  the  nonlinear  predictive  control  algorithm  for  opti¬ 
mization  is  described  as  follows.  The  task  is  to  search  for  an  optimal 
control  sequence  U{k)  {=[u{k)u{k  + 1). .  .u(k  +  M- 1)])  that  satisfies 
the  condition  that  the  objective  function  is  minimized  on  the  basis 
of  the  fuzzy  predictive  model,  the  prediction  horizon  P,  the  control 
horizon  M  (M  <  P),  and  the  regression  data  vector  x'(fc)  consisting  of 
input-output  data  at  the  current  instant  and  before: 


x  (k)  =  [u{k  —  1 ), . . . ,  u(k  —  nu ),  y{k), . . . ,  y{k  —  ny)]  (8) 


Here,  the  quadratic  objective  function  J  is  employed: 

P  M 

J  =  +  i)  -  ym(k  +  i )]2  +  Jjj  Au2(k  +  i  - 1 )  (9) 

1=1  j= 1 

where  and  Xj  are  weight  coefficient,  from  the  every  time  k,  there 
arej  control  increments  A u(/<),  ...  ,Au(k+j-  1),  y{k  +  i)  is  actual 
output,  and  ym(k  +  i)  is  model  predictive  output.  The  constraints  on 
the  system  variables  can  be  confirmed  according  to  different  system 
requirements.  The  discrete  optimization  of  the  control  sequence 
employs  a  branch-and-bound  algorithm,  which  is  an  integral  part 
of  the  programming.  The  basic  concept  of  obtaining  the  optimal 
solution  is  to  divide  the  feasible  solution  space  into  a  number  of 
small  subspaces,  and  then  to  select  the  appropriate  upper  and  lower 
boundary  functions  for  narrowing  the  solution  scope  until  the  opti¬ 
mal  solution  is  obtained  [8,21-25]. 

3.1  Determining  the  tree  structure  discrete  search  space  of  the 
control  sequence  U(k) 

The  error  e(/<)  between  the  objective  setting  value  yr(k)  of  the 
current  sampling  period  and  the  actual  system  output y(k\  and  the 
error  change  ratio  ec(/<),  may  be  calculated  as  follows: 

e(k)=yr(k)-y(k)  (10) 

ec{k)  =  e{k)  -  e(k  -  1)  (11) 

e  and  ec  are  discretized  into  E  and  EC  in  their  respective  fuzzy 
domains,  and  the  variable  du(k )  of  u(k)  is  inferred  according  to  the 
following  analytical  fuzzy  reasoning  formula: 

du(k)  =  aE  +  {  1  -a)  EC  (12) 

where  a  is  an  undetermined  coefficient. 

The  consecutive  space  of  u{k)  that  is  centered  on  u{k  -  1 )  is  given 
by 

B  =  (u(/<  -  1 )  -  abs(du),  u(k  -  1 )  +  abs(du ))  (13) 

The  consecutive  space  of  u{k)  is  then  discretized  into  the  discrete 
search  space.  Analogously,  the  discrete  search  space  of  u(k  +  i)  can 
be  obtained  according  to  yr(k  +  i)  and  y(k  +  i)  in  the  control  horizon 
(1  <i<M-  1)  based  on  each  discrete  search  point  of  u(  k  +  i  -  l).The 
control  sequence  does  not  branch  anymore  during  the  sampling 
period  beyond  the  control  horizon  (M<i<P  -  1).  Fig.  3  illustrates 
the  discrete  search  space  of  tree  structure  of  the  control  sequence 
U(k )  in  the  predictive  horizon  P  and  the  control  horizon  M.  In  this 
figure,  ut{k  +  i){t=  1, . .  .,n/<+1)  denotes  the  tth  discrete  search  point  in 
the  discrete  space  of  u  at  the  (/<  +  i)th  sample  instant.  This  algorithm 
embodies  the  branch  principle  of  the  branch-and-bound  method 

[19]- 


MPC  consists  of  the  identified  predictive  model,  the  objective 
function,  and  all  kinds  of  constraints  of  the  system  variables  [26]. 
MPC  must  find  a  control  sequence  in  the  future  limited  horizon  so 
that  the  objective  function  is  minimized,  and  the  first  control  vari¬ 
able  of  a  control  sequence  is  imposed  on  the  controlled  objective  in 
order  to  achieve  so-called  receding  horizon  optimization  control. 


3.2.  Searching  for  the  optimal  control  sequence 

The  process  of  searching  for  the  optimal  control  sequence  U{k) 
is  shown  in  Fig.  3.  In  the  discrete  search  space  of  tree  structure, 
searching  for  an  optimal  control  sequence  means  meeting  the  cri¬ 
terion  that  the  quadratic  objective  function  is  minimized.  Firstly, 
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y(k )  y(k+ 1)  y(k+M-\ )  y(k  +P-1) 


Fig.  3.  The  tree  structure  discrete  searching  space  of  control  sequence  U(k). 


(2)  While  J  is  less  than  JUpper ,  the  process  will  continue  from 
u(k  +  i  + 1)  to  u(k  +  P-  1)  on  the  basis  of  ut(/<  +  i);  J  is  assigned 
to  Jupper  (Jupper  =J )  if  J  is  still  less  than  Jupper  in  the  process. 

(3)  When  J  is  greater  than  Jupper,  the  latter  subspaces  with  respect 
to  ut{k  +  i)  are  eliminated  and  ut+i  (/<  +  /)  is  reselected. 

(4)  If  t=nki,  u(k  +  i+ 1)  must  be  reselected. 

In  this  way,  the  entire  space  is  searched  and  the  optimal  control 
sequence  U(k)  with  respect  to  the  minimum  J  is  found.  The  first- 
dimension  element  u{k)  of  U{k)  is  regarded  as  a  control  variable  for 
the  next  sampling  period  and  is  applied  to  the  controlled  objective 
so  that  a  new  sampling  output  value  y(k  + 1 )  is  obtained.  When  the 
new  input-output  data  vector  x(/<)  is  restructured,  it  is  applied  to 
receding  horizon  optimization  control  of  the  next  circulation. 


the  greedy  algorithm  is  employed  to  find  a  good  solution  in  order 
to  set  the  initial  upper  boundary  for  the  performance  index  [27]. 

It  is  then  assumed  that  u{k  +  i  - 1)  has  been  given.  The  perfor¬ 
mance  indices  of  all  discrete  control  variables  ut(k  +  i){t=  1, . .  .,n/<+1) 
are  calculated  at  the  (/<  + i)th  (i  =  1, . .  .,P-  1 )  level  according  to  Fig.  3: 

Jt(k  +  i)  =  qt\yr(k  +  i )  -  yc{k  +  i)]2  +  r,[ut(k  +  i)-u{k  +  i-  1  )]2 


(14) 


The  minimum  of  the  quadratic  objective  function  as  well  as 
u(k  +  i )  can  be  expressed  as  follows: 

Jmin(k  +  i)  =  mm  \Jt{k  +  i)}  (15) 

Jupper  is  the  initial  upper  boundary  of  performance  index  for 
searching  for  the  optimal  control  sequence  and  it  can  be  calculated 
by  the  following  formula: 

p-i 

Jupper  =  J2j(k  +  V  (16) 

z=0 

The  process  of  searching  for  the  optimal  control  sequence  U(k) 
is  shown  in  Fig.  4.  The  notation  in  the  following  steps  corresponds 
with  that  in  Fig.  4. 


(1 )  It  is  assumed  that  ut(k  +  i)  (i  =  1, . . .,  P -  1 )  is  based  on  u{k  +  i  -  1 ) 
and  J  is  calculated  as  follows: 


k+i 

J  =  5>[yr(/) 

l=k 


min{M,k+i} 

yd)]2  +  J!  r,Au(,)2 

l=k 


(17) 


(b) 


(c) 


k  +/  +1 


4.  Simulation  results 

This  section  presents  numerical  simulations  to  illustrate  the 
validation  of  the  proposed  predictive  control  developed  for  the  tem¬ 
perature  of  SOFC  stack  based  on  an  improved  T-S  fuzzy  model.  From 


Table  1 

Parameters  of  the  SOFC  stack  used  in  the  physical  model. 


Item 

Value 

Cell  number 

60 

Cell  active  area  (m2) 

0.01 

Heat  capacity  per  cell  area  (J  K_1  nrr2) 

7000 

Arm  power  (kW) 

3.5 

Anode  total  volume  (m3) 

0.005 

Cathode  total  volume  (m3) 

0.001 

Gas  inlet  flow  constant 

0.3629  x  10“5 

Anode  gas  outlet  flow  constant 

0.0544  x  10“5 

Cathode  gas  outlet  flow  constant 

0.2118  x  10“5 

Air  heat  capacity  (J  kg-1  K_1 ) 

1004 

Air  density  (kg  itt3 ) 

1.23 

Oxygen  density  (kg  m-3 ) 

1.36 

Nitrogen  density  (kg  m-3 ) 

1.19 

Air  gas  constant  (kg  nrr3 ) 

286.9 

Oxygen  gas  constant  (kg  m-3) 

259.8 

Nitrogen  gas  constant  (kg  nr3 ) 

296.8 

Hydrogen  gas  constant  (kg  nrr3 ) 

4.1243  x  103 

Oxygen  molar  mass  (kg  mol-1 ) 

0.032 

Nitrogen  molar  mass  (kg  mol-1 ) 

0.028 

Hydrogen  molar  mass  (kg  mol-1 ) 

2.016  x  10“5 

Faraday’s  constant 

96485.34 

Universal  gas  constant  (kg nr3) 

8.31451 

Fig.  4.  Demonstration  of  the  searching  process. 


Fig.  5.  Load  current  curve. 
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(a)  actual  output  and  model  predictive  output  of  stack 
temperature, and  RMSE  =  0.925 


above-mentioned  Sections  2  and  3,  the  model  predictive  control 
is  designed  based  on  the  improved  T-S  model  which  is  set  up  at 
first.  The  SOFC  physical  model  (see  Appendix  B)  described  in  Ref. 
[19]  is  used  to  obtain  the  input/output  data,  and  just  replaces  the 
true  SOFC  stack.  In  this  study,  the  training  data  is  generated  with 
the  physical  SOFC  model.  In  order  to  obtain  available  identification 
data,  the  input  signals  of  the  SOFC  physical  model  were  uniformly 
random,  mainly  is  the  current  (0-60  A).  The  parameters  of  a  3.5  kW 
SOFC  stack  were  used  in  the  simulation,  and  are  given  in  Table  1. 
The  input/output  data  was  collected  from  the  simulation,  and  then 
the  fuzzy  modeling  algorithm  was  employed  to  identify  the  T-S 
fuzzy  model.  To  validate  the  improved  T-S  fuzzy  model,  it  was  used 
to  perform  dynamic  simulation  of  the  SOFC  stack.  The  SOFC  stack 
temperature  had  to  be  kept  constant  (in  general  973.15  K).  The  ini¬ 
tialization  parameters  of  the  improved  T-S  model  were  as  follows: 
number  of  fuzzy  clustering  c  =  4,  fuzzy  parameter  m  =  3,  and  learn¬ 
ing  rate  rj  =  0.2.  Simulations  were  performed  for  all  the  schemes 
with  the  following  tuning  parameters  of  the  predictive  controller: 
prediction  horizon  P=  14,  control  horizon  M  =  3,  controlled  variable 


Fig.  7.  Response  of  model  predictive  control. 


weighting  r  =  0.2.  The  load  current  was  considered  as  a  disturbance, 
as  shown  in  Fig.  5.  Tracking  curves  and  the  deviation  between  the 
actual  output  and  the  model  predictive  output  of  the  stack  temper¬ 
ature  are  shown  in  Fig.  6.  The  deviation  between  the  actual  output 
and  the  model  predictive  control  fluctuates  within  ±2 1<  and  the 
MPC  control  strategy  displays  good  control  accuracy.  Fig.  7  shows 
that  the  SOFC  stack  is  controlled  by  the  nonlinear  predictive  con¬ 
trol  algorithm  based  on  the  improved  T-S  fuzzy  model.  MPC  adjusts 
the  operating  temperature  to  the  set  value  (973.15  K),  minimizes 
the  temperature  fluctuation,  and  provides  control  with  satisfactory 
effectiveness. 

5.  Conclusions 

Thermal  management  of  SOFC  stacks  is  very  important  because 
high  temperature  and  temperature  fluctuation  can  lead  to  seri¬ 
ous  problems.  In  this  study,  a  nonlinear  MPC  algorithm  based 
on  a  T-S  fuzzy  identification  model  has  been  proposed.  Simula¬ 
tion  results  have  shown  the  MPC  method  to  be  valid  and  that  it 
gives  good  performance  by  virtue  of  the  nonlinear  predictive  con¬ 
troller. 

It  is  clear  that  a  model  of  a  nonlinear  SOFC  stack  can  be  built  on 
the  basis  of  an  improved  T-S  fuzzy  model,  and  that  this  can  be  used 
to  predict  the  temperature  responses  online.  The  stack  temperature 
can  be  controlled  so  as  to  smoothly  maintain  the  set  value,  and  the 
simulation  results  show  the  nonlinear  predictive  controller  to  be 
superior  for  this  purpose. 
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Appendix  A.  Nomenclature 


A  cell  active  area  (m2) 

Cp  specific  heat  capacity  of  gas  species  (J  mol-1  K-1 ) 

Cp  area  specific  heat  capacity  (J  K-1  m-2 ) 

F  Faraday’s  constant  (96,485  C  mol-1 ) 

H  enthalpy  (J) 

H  enthalpy  flow  (W) 

h  specific  enthalpy  (J  mol-1 ) 

/  current  (A) 

i  current  density  (A  nrr 2 ) 

i0  exchange  current  density  (A  m-2 ) 

ir  reaction  current  density  (A  m-2 ) 

m  mass  (kg) 

N  molar  flux  (mol  nrr2  s-1) 

n  number  of  moles  (mol) 

h  molar  flow  (mol  s_1) 

rii  molar  number  of  species  i  (mol) 

hc  hydrogen  combustion  molar  flow  (mol  s-1 ) 

P  pressure  (bar) 

R  universal  gas  constant  (8.314  J  mol-2  K-2) 

T  temperature  (K) 

T0  ambient  temperature  (I<) 

t  time(s) 

U  voltage  (V) 

V  volume  (m2) 

x  input  vector 

Xj  molar  fraction  of  species 

y  output  vector 
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Greek  letters 

a  symmetry  factor,  learning  factor 

P  a  T-S  factor 

v  stoichiometric  coefficient 

r]  overpotential  (V) 

Subscripts 

i  gas  species  i 

r  reaction 

Superscripts 
in  fuel  cell  inlet 

out  fuel  cell  outlet 

Appendix  B.  A  physical  model  of  an  SOFC  stack 


Inserting  Eq.  (B.4)  into  Eq.  (B.l)  and  assuming  constant  pressure, 
the  molar  balance  becomes: 


dXj  _  -  inf yin 

RTdF“n 


(B.5) 


B.2.  Energy  balance  model 

Here  only  the  enthalpy  balance  is  considered.  The  main  sources 
and  sinks  of  heat  for  an  SOFC  are  entering  and  exiting  flow,  the  heat 
generated  by  cell  reactions,  and  the  heat  lost  to  the  environment. 
An  enthalpy  balance  equation  yields: 

rll-t 

—  =  Hin  -  Hout  -  irVA  -  Hl05S  (B.6) 


In  this  SOFC  physical  model,  some  simplifications  and  assump¬ 
tions  are  made,  since  high  accuracy  is  not  necessary  for  a  physical 
model  that  serves  for  a  subsequent  control  strategy.  Any  deviations 
between  the  model  and  the  real  fuel  cell  can  be  managed  by  a 
feedback  loop  in  the  control  system  [28].  This  section  has  many 
symbols  which  can  be  defined  in  Appendix  A.  The  proposed  SOFC 
stack  physical  model  has  the  following  assumptions: 

(1 )  Stack  is  fed  with  hydrogen  and  air;  the  fuel  processor  dynamics 
are  not  included. 

(2)  A  uniform  gas  distribution  among  the  cells  is  assumed,  since 
there  is  a  small  deviation  of  the  gas  distribution. 

(3)  There  is  no  heat  transfer  among  the  cells.  Each  cell  has  the  same 
temperature  and  current  density. 

(4)  There  is  no  heat  exchange  between  the  stack  and  the  ambient 
environment. 

(5)  The  channels  that  transport  gases  along  the  electrodes  have  a 
fixed  volume  and  small  dimension,  so  there  is  a  constant  pres¬ 
sure  in  the  stack. 


The  term  irVA  represents  the  electrical  power  generated.  The 
relationship  between  enthalpy  and  temperature  can  be  expressed 
as: 


dH  A.  dT 
dt  p  dt 


(B.7) 


The  entering  enthalpy  flow  also  contains  the  heat  generated  by 
the  combustion  of  hydrogen,  hence  the  total  entering  enthalpy  flow 
is  equal  to  the  enthalpy  associated  with  the  ambient  air  plus  the 
entering  hydrogen  flow: 


Hin  =  naiAir(7'o)  +  (nH2  +nfj2)/iH2(rH2)  (B.8) 

For  simplification  Th2  =  To  is  assumed  without  causing  any 
meaningful  deviation. 

The  total  outlet  molar  flow  will  be  larger  than  or  equal  to  the 
entering  one  and  the  expression  of  this  is 

Hout  =  'Y)h°uthi(T)  (B.9) 


B.l.  Mass  balance  model 


For  a  generic  species  i,  the  dynamic  mole  balance  is 
^  =  ni',x'n  -  h0UtXi  +  v~ 


(B.l) 


The  stoichiometric  factor  vt  indicates  how  many  moles  of  the 
species  are  produced  or  consumed  for  each  mole  of  electrons 
transferred.  The  anodic  and  cathodic  reactions  in  the  SOFC  are, 
respectively: 


H2  +  02-^  H20  +  2e“ 


and  the  outlet  molar  flow  for  species  i ,  h?ut  can  also  be  expressed 
as: 

„r  =  n;«  +  vl(^2  +  ^)+x,|^f  (B.10) 

where  h|.n  is  the  entering  molar  flow,  v,- (h^2  +  ( irA/2F ))  is  the  molar 
flow  associated  with  hydrogen  reaction  and  combustion,  and  the 
term  xz  ( pV/RT 2)  (dT/dt)  is  related  to  gas  thermal  expansion.  For  the 
sake  of  simplicity,  it  is  assumed  that  Hloss  =  0  [28]. 

B.3.  Temperature  control  dynamic  model 


l/202  +  2e_  02~  (B.2) 

then  vq2  =  -1/2,  vh2  =  -1,  vH2o  =  1,  and  vn2  =  0. 

According  to  the  Butler-Volmer  equation  [29]: 

ir  =  i0(e^nF/RT)rj  _  e-(l -a^nF/RT)^  (B.3) 


When  pressure  is  constant,  the  total  outlet  molar  flow  depends 
only  on  the  transient  reaction  rate  and  temperature  and  can  be 
calculated  using  the  following  relation: 


EirA  pV  dT 

v‘T  +  Rf2W 

i 


(B.4) 


The  term  ^\vj( irA/F )  represents  the  algebraic  sum  of  gas  reac¬ 
tion  molar  flow,  and  ( pV/RT 2)  (dT/dt)  reflects  thermal  expansion. 


Having  assumed  that  temperature  is  uniform  in  a  stack,  Eq.  (B.6) 
determines  its  dynamics.  According  to  simultaneous  Eqs.  (B.6)  and 
(B.7),  it  can  be  obtained  that 

AcJL  =  Hin  -  Hout  -  irVA 
dt 

=  (H/n(T)  -  H0Ut(T))  -  (Hin(T)  -  Hin(To))  -  irVA  (B.ll) 


The  terms  (Hin(T)  -  H0Ut(T))  and  (H/n(T)  -  Hin(T0)),  respectively, 
represent  the  reaction  heat  and  the  sensible  heat. 

In  Eq.  (B.ll)  it  assumed  that  the  specific  heat  capacity  for  all 
the  species  is  approximately  the  same  on  condition  that  the  pres¬ 
sure  is  constant.  This  is  adequate  for  the  control-oriented  modeling. 
According  to  Eq.  (B.ll),  the  temperature  control  model  can  be  built 
with  MATLAB. 
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B.4.  Temperature  control  model  realization  in  MATLAB 

This  physical  MATLAB  model  has  two  input  variables,  air  flow 
and  H2  flow,  and  one  output  variable,  T.  A  virtual  3.5  kW  SOFC  stack 
is  used  in  the  simulation.  The  stack  consists  of  60  cells  with  anode 
and  cathode  gases  in  cross-flow  and  cell  active  area  0.01  m2.  The 
physical  model  replaces  the  real  SOFC  stack  to  generate  the  simula¬ 
tion  data  required  for  the  modified  T-S  fuzzy  model.  The  enthalpy 
expressions  can  be  acquired  from  the  published  literature  [30]  and 
the  relationship  between  enthalpy  and  temperature  for  different 
gases  is  reflected  in  these  expressions.  The  Enthalpy _out  block  and 
Enthalpy _in  block  can  be  expressed  as  follows  (w  expresses  molar 
flow): 

For  Enthalpy _in  block: 

hJn  =  w_air  x  h_air  +  w_H2  x  h_H2  +  w_V  x  h_V  (B.12) 

where 

h_air  =  -1.0947  e4  +  32.50T_air 


h_H2  =  -0.9959  e4  +  30.73T_fuel 

hM  =  -25.790  e4  +  42.47T_fuel  (B.13) 

w_FI2  =  w_fuel  x/Lmoelra  x  2.016/(/Lmoelra  x  2.016  +  18.02) 

w_V  =  w_fuel  -  w_H2  (B.14) 

For  the  Enthalpy _out  block: 

h_out  =  w_02  x  h_02  +  w_N2  x  h_N2  +  w_H2  x  h_H2  +  w_V  x  h_V 

(B.15) 

and 

h. 02  =  -1.2290e4  +  35.12T 


h_N2  =  -1 .0590  e4  +  31 .401 


h_FI2  =  —0.9959  e4  +  30.73T 


h_V  =  -25.790  e4  +  42.47T  (B.16) 

Based  on  the  physical  model,  it  can  be  established  that  the  load 
signal  change  can  vary  the  SOFC  stack  temperature. 
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